function [cost_measure spikeout mfr all_target_mfr] = spnet( conn_data_in, optional_pattern, all_target_mfr)

if nargin > 0
    if size(conn_data_in,2)==1
        conn_data_in = reshape(conn_data_in, 6, 3);
    end
end
start_save_sec = 0;
stop_sec = 1;


if nargin > 1
    r=optional_pattern;
else
     r=[.122 .33 .12 .47 1 .42 .03 .43 .31 .31 .45; ...
     0.2848    0.1522    0.1166    0.2534    0.4529    0.4443    0.1073    0.1756    0.2698    0.3289    0.7500; ...
     0.1543    0.2126    0.5000    0.1736    0.1401    0.0113    0.0370    0.2063    0.2159    0.2477    0.1655; ...
     0.1009    0.1317    0.1265    0.0578    0.0548    0.2500    0.1192    0.0408    0.0170    0.0298    0.1115];
%r=[    0.0321    0.4475    0.2013    0.9687    0.1278    0.2822    0.2924    0.1789    0.0484    0.2318    0.7734;...
%    0.7965    0.8375    0.8627    0.0655    0.6800    0.4220    0.9668    0.1903    0.9100    0.5607    0.0139;...
%    0.0341    0.7960    0.2558    0.2295    0.0283    0.6533    0.7794    0.0992    0.1220    0.5733    0.2672;...
%    0.0185    0.9369    0.0597    0.0331    0.8519    0.9216    0.0863    0.9099    0.3537    0.3028    0.9321];

    all_target_mfr = [...
         0.0007    0.0091   14.0048   37.5236   47.7772   37.5828   13.5727    0.0255    0.0067    0.0053    0.0084;...
         29.3920    7.6146    0.0011    0.0028    0.0073    0.0077    0.0003    0.0069   18.3303   35.9531   42.0335;...
         18.9614   25.8034   24.2515   11.4287    0.0012         0         0    0.0032    0.0051    0.0224    5.1803;...
           0.0009    0.0068    0.2328    4.4247    7.4911    8.7554    5.4898    0.3377         0         0    0.0006 ...
            ];


        
    epoch_ms = 1000;
    npatterns = 4;
    isi_ms = floor(epoch_ms/npatterns);

        
end

search_params = [1 3 4];
conn_data = zeros(6,9);
    

conn_arg = 1;

%function spikeout = spnet( stop_sec, start_save_sec,  conn_data)
% spikeout will be all spikes between start_savfe_sec and stop_sec.

global s sd syntype D N syndata v I Ne IN A AI vhist ehist ampahist nmdahist gabaahist gababhist Achand sum_tempi sum_tempe

draw = true;
learning = false;
  

spikeout=[];

% spnet.m: Spiking network with axonal conduction delays and STDP
% Created by Eugene M.Izhikevich.                February 3, 2004
% Modified to allow arbitrary delay distributions.  April 16,2008
% Rewritten by Jeffrey L. McKinstry, November, 2009 to allow anatomy.

D=5;                  % maximal conduction delay 
% excitatory neurons   % inhibitory neurons      % total number 
Ne=22;                Ni=22;                   N=Ne+Ni;
% Copyright 2009 Neurosciences Research Foundation, Incorporated

% set up group sizes.
IN = 1:11;
A = max(IN)+1:Ne;
Achand = max( A ) + 1:floor(Ne+Ni/2);
AI = max(Achand)+1:N;
N = max(AI);


% neuron parametes:
a=[0.02*ones(Ne,1);    0.1*ones(Ni,1)];
d=[   8*ones(Ne,1);    2*ones(Ni,1)];

nmda_ampa_ratio=ones(N,1);
nmda_ampa_ratio(Achand)=.2;
nmda_ampa_ratio(AI)=.12;

gabab_gabaa_ratio=0.0*ones(N,1);
gabab_gabaa_ratio(Achand)=0;
gabab_gabaa_ratio(AI)=0;

% Take special care not to have multiple connections between neurons
s = cell(D,1);
sd = cell(D,1);
syntype = cell(D,1);
for i = 1:D
    s{i} = sparse(N,N);
    syntype{i} = sparse(N,N);
    sd{i} = sparse(N,N);
end



% connections
%r=rand(1,12)*.5;  r(5)=1

patterns{1} = r(1,:);
patterns{2} = r(2,:);
patterns{3} = r(3,:);
patterns{4} = r(4,:); %[0 0 0 0 0 0 0 0 0 0 0 0];



%patterns{1} = [0 0 0 .5 1 .5 0 0 0 0 0 0];
%patterns{2} = [0 0 0 .5 1 .5 0 0 0 0 0 0];
%patterns{3} = [0 0 0 0 0 0 0 0 0 0 0 0];
%patterns{4} = [0 0 0 0 0 0 0 0 0 0 0 0];

%patterns{1} = [1 1 1 1 1 1 0 0 0 0 0 0];
%patterns{2} = [0 0 0 1 1 1 0 0 0 1 1 1];
%patterns{3} = [1 1 1 0 0 0 1 1 1 0 0 0];
%patterns{4} = [0 0 0 0 0 0 1 1 1 1 1 1];


syndata = [];


i=1;


radius=.4;
conn_data(i,:)=[0.5*radius 0 1 radius 3 3 .05 0 10]; 
if nargin >= conn_arg,
    % merge the supplied parameters with defaults.
    conn_data(i,search_params)=conn_data_in(i,:);    
end
conn = conn_data(i,:);i=i+1; 
connect(IN,A,'gauss', 'std', conn(1), 'initoffsetw', conn(2), 'initrangew', conn(3), 'radius', conn(4), 'delaymin', conn(5), 'delayrange', conn(6), 'lrate', conn(7), 'minw', conn(8), 'maxrangew', conn(9));

% This might be helpful in combination with short-term syn. plasticity
% That reduces this value on inhib neurons quickly.

% radius=.4;
% conn_data(i,:)=[0.5*radius 0 1 radius 3 3 .05 0 15];
% if nargin >= conn_arg,
%     % merge the supplied parameters with defaults.
%     conn_data(i,search_params)=conn_data_in(i,:);    
% end
% conn = conn_data(i,:);i=i+1; 
% connect(IN,AI,'gauss', 'std', conn(1), 'initoffsetw', conn(2), 'initrangew', conn(3), 'radius', conn(4), 'delaymin', conn(5), 'delayrange', conn(6), 'lrate', conn(7), 'minw', conn(8), 'maxrangew', conn(9));

radius=0.4;
conn_data(i,:)=[0.5*radius 0 1 radius 1 2 .05 0 20]; 
if nargin >= conn_arg,
    % merge the supplied parameters with defaults.
    conn_data(i,search_params)=conn_data_in(i,:);    
end
conn = conn_data(i,:);i=i+1; 
connect(A,AI,'gauss', 'std', conn(1), 'initoffsetw', conn(2), 'initrangew', conn(3), 'radius', conn(4), 'delaymin', conn(5), 'delayrange', conn(6), 'lrate', conn(7), 'minw', conn(8), 'maxrangew', conn(9));


radius=0.4;
conn_data(i,:)=[radius 0 1 radius 1 2 .05 0 20]; 
if nargin >= conn_arg,
    % merge the supplied parameters with defaults.
    conn_data(i,search_params)=conn_data_in(i,:);    
end
conn = conn_data(i,:);i=i+1; 
connect(A,Achand,'gauss', 'std', conn(1), 'initoffsetw', conn(2), 'initrangew', conn(3), 'radius', conn(4), 'delaymin', conn(5), 'delayrange', conn(6), 'lrate', conn(7), 'minw', conn(8), 'maxrangew', conn(9));

% make chandelier cells excitatory.
radius = length(A)/2;
not_used = pi;
conn_data(i,:)=[not_used 0  4 radius 1 2 0.0 0 10];
if nargin >= conn_arg,
    % merge the supplied parameters with defaults.
    conn_data(i,search_params)=conn_data_in(i,:);    
end
conn = conn_data(i,:);i=i+1; 
connect(Achand,A,'cos', 'std', conn(1), 'initoffsetw', conn(2), 'initrangew', conn(3), 'radius', conn(4), 'delaymin', conn(5), 'delayrange', conn(6), 'lrate', conn(7), 'minw', conn(8), 'maxrangew', conn(9));

radius = length(A)/2;
conn_data(i,:)=[not_used 0 4 radius 1 1 0.0 0 20]; 
if nargin >= conn_arg,
    % merge the supplied parameters with defaults.
    conn_data(i,search_params)=conn_data_in(i,:);    
end
conn = conn_data(i,:);i=i+1; 
connect(AI,A,'cossurround', 'std', conn(1), 'initoffsetw', conn(2), 'initrangew', conn(3), 'radius', conn(4), 'delaymin', conn(5), 'delayrange', conn(6), 'lrate', conn(7), 'minw', conn(8), 'maxrangew', conn(9));


mfr=zeros(N,1);




%s=[6*ones(Ne,M);-5*ones(Ni,M)];         % synaptic weights

ltp_decay = exp(-1/20);
ltpcurve = .1*cumprod([ 1 ltp_decay*ones(1,49)]);
ltd_decay = exp(-1/20);
ltdcurve = -.12*cumprod([ 1 ltd_decay*ones(1,49)]);
stdpcurve = [ fliplr(ltpcurve)  ltdcurve 0];
  
if draw
    mv = zeros(N,1); % for visualization
    vhist = zeros(N,1000);
    ihist = zeros(N,1000);
    ehist = zeros(N,1000);
    ampahist = zeros(N,1000);
    nmdahist = zeros(N,1000);
    gabaahist = zeros(N,1000);
    gababhist = zeros(N,1000);
    itothist = zeros(N,1000);    
end


% restart for each simulation.
%rand('seed',1)
tau_ampa  = exp(-1/5);
tau_nmda  = exp(-1/150);
tau_gabaa = exp(-1/6);
tau_gabab = exp(-1/150);


v = 0.11*rand(N,1);                      % initial values modified for mfr model

MAXSPIKESPERSEC= N/2*1000;
firings=zeros(MAXSPIKESPERSEC,2);

firehist=cell(N,1);

I=zeros(N,1);
Igabaa = I;Igabab = I; Iampa = I;Inmda = I;
Ie=zeros(N,D);
Ii=zeros(N,D);

numspikes = 0;

snew = [];
for i=1:D
    snew = [snew s{i}];
end
snew=full(snew);
se=snew(1:max(Achand),:);
si=snew(max(Achand)+1:N,:);


for sec=1:stop_sec                      % simulation of 1 day
  for t=1:1000                          % simulation of 1 sec
    
    %randindex = ceil(N*rand);
    %I(randindex,1)=I(randindex,1)+20;                 % random thalamic input 
%    I(IN,1) = I(IN,1) + 20*patterns{ 1+floor((t-1)/250)};
    
    efired = sqrt(max(0.1, v' ));

    ifired = efired(max(Achand)+1:N);
    efired = efired(1:max(Achand));

    
    % deliver all the synapses to post targets with correct delays.
    % I is a shift register for each neuron.
    %for i = 1:D
     %  Ie(:,i) = Ie(:,i) + (spfired(1:max(Achand))*s{i}(1:max(Achand),:))';
     %  Ii(:,i) = Ii(:,i) + (spfired(max(Achand)+1:N)*s{i}(max(Achand)+1:N,:))';
    %end
     %Ie = Ie + reshape( (spfired(1:max(Achand))*snew(1:max(Achand),:)) , N,D);
     %Ii = Ii + reshape( (spfired(max(Achand)+1:N)*snew(max(Achand)+1:N,:)) ,N,D);
     Ie = Ie + reshape( efired*se , N,D);
     Ii = Ii + reshape( ifired*si ,N,D);


    Igabaa = Igabaa*tau_gabaa + Ii(:,1);
    Igabab = Igabab*tau_gabab + gabab_gabaa_ratio.*Ii(:,1);

    Iampa = Iampa*tau_ampa + Ie(:,1);
    Inmda = Inmda*tau_nmda + nmda_ampa_ratio.*Ie(:,1);

    
    xx = (v)/60;      % modified for mfr model
%    xx = (-80-v)/60;
	xx = xx.*xx;
	NMDAgate = xx./(1+xx);
    NMDAgate(v<0)=0;  % modified for mfr model
    
    I = Igabaa + Igabab + Iampa;% + NMDAgate.*Inmda;
    I(IN) = I(IN) + 10*patterns{ 1+floor((t-1)/isi_ms)}';

    alpha = exp(-1/10);
    v = alpha*v + (1-alpha)*I;

    %I= [I(:,2:D) zeros(N,1)];  % rotate.     
    Ie= [Ie(:,2:D) zeros(N,1)];  % rotate.     
    Ii= [Ii(:,2:D) zeros(N,1)];  % rotate.     
    
    
  
  
    if draw
        vhist(:,t)=v;
        
        tau = exp(-1/200);
        mfr=tau*mfr + (1-tau)*v;
        mfr=mfr*tau;

        mv = tau*mv + (1-tau)*v;
        if mod(t,20)==0
            figure(3);
            plot([v(A) I(A) I(IN)])
            hold on
            drawnow;
            
            %figure(4)
            %plot(mv);
            %drawnow;
        end
    end
    
  end;
  
  
  disp(['t=' num2str(sec)])
  
  
  
  
end;
hold off

% Calculate cost measure

% build the desired mfr profile.


% all_target_mfr = zeros(npatterns,length(A));
% for pat = 1:npatterns
% 
%     std=length(A)/10;
%     radius=floor(length(A)/2);
%     target_mfr=(gaussian_1d(std,-radius:radius));
%     target_mfr=target_mfr*(1/max(target_mfr)).*30; % winner fires at 30 hz.
%     [m max_index]=max( patterns{pat} );  % shift s.t. winner is in the middle.
%     target_mfr=target_mfr(1:length(A));
%     % assumes length(A) = length(r);
%     mid = ceil(length(A)/2);
%     shift = abs(max_index - mid)
%     if max_index > mid
%         % shift right
%         target_mfr=[ target_mfr(end-shift+1:end) target_mfr(1:end-shift)];
%     else
%         % mid >= max_index
%         target_mfr=[ target_mfr(shift+1:end) target_mfr(1:shift)];    
%     end
% 
%     all_target_mfr(pat,:)=target_mfr;
% end
    
% bin the spike data by pattern.
if ~isempty(spikeout)    
    S=sparse( spikeout(:,2), spikeout(:,1), 1, N , stop_sec*1000 );

    % Gather mean rates for stats.

    t=start_save_sec*1000 + 1;
    count = zeros(1,npatterns);
    mfr = zeros(npatterns,N);
    for i = start_save_sec*1000:epoch_ms:stop_sec*1000-1
        for p = 1:npatterns
            mfr(p, : ) = mfr(p, :) + (1000/isi_ms)*full(sum( S( :, t+20:(t+isi_ms-1)), 2))';
            count(p)=count(p)+1;
            t=t+isi_ms;
        end
    end
    mfr=mfr./(repmat( count', 1, size(mfr,2) ) + .000001);


    % Look at sparsness measures
    %     
    %     j= 1;
    %     for i = max(IN)+1:Ne
    %         lifetime_sparseness(j) = sparseness( mfr(:,i), 40 ); % 40 Hz is max
    %         j= j + 1;
    %     end
    %     disp(['Mean lifetime sparseness: ' num2str( mean(lifetime_sparseness))]);
    %     
    %     for i = 1:patterns
    %         population_sparseness(i) = sparseness(mfr(i,max(IN)+1:Ne), 40);
    %     end
    %     disp(['Mean population sparseness: ' num2str( mean(population_sparseness))]);

        % save the best
    %    d= sqrt( (mean(population_sparseness)-0.5)^2 )
    %d= sqrt( (mean(population_sparseness)-0.5)^2 + (mean(lifetime_sparseness)-0.5)^2)

    cost_measure = sqrt( sum(sum((mfr(:,max(IN)+1:Ne) - all_target_mfr ).^2 )))

else
    cost_measure = 1e10;
end


stop


%
% possible arguments:
%  arbor_type == 'gauss'
%  'std', 'minw', 'maxrangew', 'dist', 'delayrange', 'lrate';
%
function connect( from_indices, to_indices, arbor_type, varargin)

global s sd syntype N D syndata


L = length( syndata);


for k= 1 :2: size(varargin,2) 

        switch lower( varargin{k} )
          case {'std'}
            std = varargin{k+1};
          case 'initoffsetw'
            minw = varargin{k+1};
          case 'initrangew'
            maxw = varargin{k+1};
          case 'radius'
            radius = round(varargin{k+1});
          case 'delaymin'
            delaymin = round(varargin{k+1});
          case 'delayrange'
            delayrange = round(varargin{k+1});
            if delaymin + delayrange-1 > D
                disp('delay > max delay in connect. Using max.');
                delayrange
                delaymin
                delayrange = D-delaymin + 1;
            end
          case 'lrate'
            syndata(L+1).lrate = varargin{k+1};
          case 'minw'
            syndata(L+1).min = varargin{k+1};
          case 'maxrangew'
            syndata(L+1).max = syndata(L+1).min + varargin{k+1};
          otherwise
            disp(['Unknown param.' varargin{k}])
        end

end

% Assume wraparound.

% clip radius to not touch target neurons more than once.
if radius > length(to_indices)/2-.5, radius = floor(length(to_indices)/2-.5);end

p=[to_indices to_indices to_indices]; %randperm(N);
% create the list of indices for each neuron; project topographically.
step = length(to_indices)/length(from_indices);
rel_index = .5;

if strcmp(arbor_type, 'gauss')
    w = minw+maxw *std*sqrt(2*pi)* gaussian_1d(std, -round(radius):round(radius) );
elseif strcmp(arbor_type, 'surround')
    % surround
    x=-round(radius):round(radius);
    w=std*sqrt(2*pi)* abs(x).*gaussian_1d(std, x );
    w = minw+maxw*(1/max(w))*w
elseif strcmp(arbor_type, 'cos')
    w=maxw*cos(-1.5*pi:1.5*pi/radius:1.5*pi)
    w(w<0)=0;
elseif strcmp(arbor_type, 'cossurround')
    w=maxw*cos(-1.5*pi:1.5*pi/radius:1.5*pi)
    w(w<0)=w(w<0)*(sum(w(w>0))/abs(sum(w(w<0))));
    w(w>0)=0;
end
for i=from_indices
    my_center = round(rel_index*step);
    rel_index=rel_index + 1;
    myoffset = length(to_indices) + my_center;

    target_indices = p( myoffset + (-round(radius):round(radius)) );
    
    delays = delaymin + floor( delayrange*(abs(-round(radius):round(radius))/(round(radius)+.001)));

    for j = 1:length(target_indices)
        s{delays(j)}(i,target_indices(j)) = w(j); % + noise.
        syntype{delays(j)}(i,target_indices(j)) = length( syndata );
        sd{delays(j)}(i,target_indices(j)) = 0.0000001; % + noise.        
    end

    
end;





